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Capon’s  maximum  likelihood  (ML)  method  has  been  used  with  some  success  in  matched  field 
processing  for  source  range  and  depth  estimation.  One  reason  for  the  interest  in  the  ML  is  that 
it  is  data  adaptive;  that  is,  it  adapts  to  the  actual  noise  field  present  rather  than  requiring  an 
a  priori  estimate  of  the  noise  component  for  prewhitening.  When  modal  noise  is  present  the 
ML  can  become  sensitive  to  any  deviations  from  the  unperturbed  case  (i.e.,  from  the  model)  as 
would  be  introduced  by  phase  errors  or  model  parameter  errors.  Using  a  dimensionality- 
reduction  procedure  a  more  stable  data  adaptive  method,  the  “reduced”  ML  (RML),  is 
obtained.  The  RML  is  compared  here  with  the  ML  on  simulated  data  from  a  21-sensor  array 
in  a  Pekeris  waveguide  supporting  eight  normal  modes.  Under  modal  noise  conditions  the 
RML  provides  a  significant  improvement  over  ML  when  phase  errors  occur.  Although  the 
deviation  from  the  model  conside^d  here  is  that  caused  by  phase  errors,  the  nature  of  the 
perturbation  is  not  important  sii  .he  sensitivity  of  ML  is  not  to  any  special  type  of 
perturbation. 
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INTRODUCTION 

A  number  of  recent  articles  have  considered  the  topic  of 
matched  field  processing,  in  which  the  plane  wave  steering 
vector  used  in  bearing  estimation  is  replaced  by  a  more  com¬ 
plicated  vector  that  models  what  appears  at  the  array  when 
there  is  a  source  at  a  particular  range  and  depth  in  a  wave¬ 
guide.' Most  published  work  to  date  deals  with  a  range- 
independent  propagation  model  and  employs  a  normal 
mode  representation  of  the  pressure  field.  When  a  conven¬ 
tional  matched  filter  approach  to  range-depth  estimation  is 
taken,  the  resulting  ambiguity  surface  often  contains  several 
false  peaks.  To  obtain  smoother  surfaces  one  can  use  various 
"mode  space"  beamforming  methods4  ';  recently,  the  use  of 
Capon’s  maximum  likelihood  ( ML)  method6  7  has  also  been 
considered. 

It  has  been  shown  for  the  plane  waves  case"  that  Capon’s 
ML  can  become  quite  sensitive  to  deviations  from  the  unper¬ 
turbed  model  case  in  the  presence  of  a  correlated  noise  field, 
such  as  spherical  isotropic  noise  with  an  oversampled  array. 
For  the  more  complicated  case  of  normal  mode  propagation 
in  a  waveguide,  the  ML  again  exhibits  sensitivity  to  pertur¬ 
bations  in  the  data,  particularly  in  the  presence  of  modal 
noise.  “Sector-focused  stability”  (SFS)  was  introduced  in 
Ref.  8  to  stabilize  the  ML  in  the  plane  wave  case.  The  SFS 
involves  the  use  of  a  square  matrix  whose  size  can  be  consid¬ 
erably  less  than  the  number  of  sensors,  and  therefore  reduces 
computation  cost  and  numerical  errors  while  improving  sta¬ 
bility.  When  the  philosophy  behind  SFS  is  applied  to  the  case 


of  normal  mode  propagation  in  a  waveguide  we  obtain  the 
“reduced”  ML  (RML). 

To  understand  why  the  ML  is  unstable  in  the  cases  con¬ 
sidered  we  examine  the  eigenvector/eigenvalue  decomposi¬ 
tion  of  the  cross  spectral  matrix  obtained  from  the  data  and 
expand  the  ML  estimator  as  the  reciprocal  of  a  sum  of  terms, 
one  for  each  eigenvector.  Most  of  the  terms  in  the  denomina¬ 
tor  of  the  ML  estimator  will  be  zero  at  the  correct  range  and 
depth,  leading  to  a  large  peak  in  the  ML  estimator.  When 
there  are  no  phase  errors,  several  terms  in  the  sum  vanish 
identically  across  the  ambiguity  surface.  When  phase  errors 
are  present  these  terms  are  slightly  bigger  than  zero  every¬ 
where  and  the  value  of  the  denominator  at  the  true  range  and 
depth  increases,  leading  to  a  drop  in  peak  height.  By  remov¬ 
ing  the  offending  small,  but  positive,  terms  one  can  restore 
the  peak  and  stabilize  the  ML  estimator;  when  so  modified 
the  estimator  becomes  the  RML.  It  is  not  necessary  to  calcu¬ 
late  each  eigenvector  and  each  term  separately;  a  dimension¬ 
ality-reduction  procedure  accomplishes  this  and  avoids  the 
inversion  of  the  full-sized  matrix  as  well. 

The  paper  is  organized  as  follows:  Sec.  I  formulates  the 
problem  in  mathematical  language,  presenting  the  normal 
mode  model  for  the  vector  of  pressure  field  samples,  intro¬ 
ducing  the  theoretical  cross  spectral  matrix  to  be  used 
throughout,  and  touching  briefly  on  two  linear  estimators, 
the  conventional  matched  filter  method  (CONV)  and  the 
mode  space  version  of  the  conventional  estimator,  called 
here  the  modified  conventional  method  (MCONV).  In  Sec. 
II  we  present  Capon’s  ML  method,  discuss  proper  normali- 
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nation  procedures,  introduce  the  eigenvector/eigenvalue  de¬ 
composition.  and  consider  the  source  of  the  instability.  In 
Sec.  Ill  we  develop  the  RML  and  present  alternate  methods 
for  its  calculation.  Section  IV  contains  a  discusion  of  our 
simulations  and  results.  Section  V  contains  comments  and 
conclusions.  We  include  three  appendices.  The  first  extends 
our  discussion  of  the  sensitivity  problem  to  local  normal 
mode  models  for  the  range-dependent  case,  and  the  second  is 
concerned  with  similar  sensitivity  in  the  use  of  the  parabolic 
approximation.  The  third  appendix  contains  a  mathematical 
derivation  relevant  to  the  normalization  of  the  ML  estima¬ 
tor. 

I.  FORMULATION  AND  BACKGROUND 

We  consider  a  vertical  array  of  .V  sensors  in  an  acoustic 
waveguide  supporting  M  <  .V  discrete  propagating  modes. 
We  shall  be  concerned  here  with  the  case  studied  by  Buck¬ 
ingham/'  the  so-called  "low-loss"  situation,  in  w'hich  the 
main  contribution  to  the  ambient  noise  field  is  discrete  mod¬ 
al  propagation  over  long  distances  of  acoustic  energy  due  to 
independent,  randomly  distributed,  surface  sources  (see 
also  Kuperman  and  Ingenito1").  For  the  case  we  consider 
here  the  single-snapshot,  narrow-band,  data  vector  x  (A  by 
1 )  has  the  algebraic  form 

x  =  a()Us  +  Uy  +  -q  .  (1) 

Flere,  U  is  the  N  by  M  matrix  whose  («,m)th  entry  is 
U„,„  =  U,„  (z„ ).  which  is  the  evaluation,  at  the  depth  z„  of 
the  nth  sensor,  of  the  mth  mode  amplitude  function  (see 
Kuperman  and  Ingenito,"1  p.  1990).  The  entries  of  s  and 
(random  vector)  y  represent  the  excitation  of  the  various 
mode  amplitude  functions  by  the  signal  and  the  aggregate 
noise  sources,  respectively.  Also,  a„  is  a  complex  Gaussian 
random  amplitude  and  q  is  a  random  white  (nonmodal) 
noise  vector. 

For  the  case  of  range-independent  sound  speed  and  den¬ 
sity,  and  assuming  cylindrical  symmetry,  the  pressure  field 
at  depth  z  due  to  a  single  source  at  depth  z,  and  range  r  is 
given  by 

M 

p(r,zt,z)  =  £  U,„(z)s,„(r,zJ  ,  (2) 

tn  -  1 

where  the  U,„(z)  are  the  samples,  at  the  sensor  depth  z,  of 
the  m th  modal  eigenfunction  that  arises  as  a  solution  to  the 
depth  variable  Sturm-Liouville  boundary  value  problem," 
and  the  sm{r,zs )  are  the  modal  amplitudes  given  by 

sm{r,zj  =  Um(zjexp(ikmr)(2ir/kmr)'n- .  (3) 

The  matrix  U  has  entries  Unm  =  Um  (z„ ),  wherez,„..,zv  are 
the  depths  of  the  N  sensors  and  the  vector  s  has  the  entries 
s,„(r,zs).  Here,  k2m  is  the  eigenvalue  associated  with  the  mth 
mode.  The  noise  vector  Uy  represents  modal  noise  created 
by  the  excitation  of  the  modal  structure  by  distributed  inde¬ 
pendent  sources  of  ambient  acoustic  energy,  such  as  wave 
action  on  the  surface  and/or  many  range-separated  surface 
ships  (see  Kuperman  and  Ingenito,1"  p.  1990). 

When  the  sound  speed  or  the  density  is  range  depen¬ 
dent,  one  can  modify  the  normal-mode  analysis  to  obtain  a 
theory  based  on  local  normal-mode  structure,"  12  or  use  a 
parabolic  approximation. 1 1 1n  either  case  the  model  given  by 

2494  J  Acoust  Soc  Am.,  Vol  87,  No.  6.  June  1990 


Eq.  ( 1 )  applies.  We  consider  these  two  approaches  to  the 
range-dependent  case  in  Appendices  A  and  B. 

For  x  given  by  Eq.  (1),  the  cross-spectral  matrix 
(CSM)  is/?  =  (xx"), with//  denotingconjugate transposi¬ 
tion  and  (  )  time  averaging;  R  is  the  average,  over  the  avail¬ 
able  snapshots,  of  the  dyad  matrices  xx",  and  is  given  (ap¬ 
proximately)  by 

R  =  \ar(a(l)Usns"U"  +  UQU"  +p2I ,  (4) 

where  Q  is  the  CSM  associated  with  the  modal  noise  field 
( i.e„  Q  =  (yy11 ) )  and  p 2  is  the  white  noise  power.  The  ma¬ 
trix  Q  may  or  may  not  be  a  diagonal  matrix.  This  depends  on 
the  geometric  distribution  of  random  sources  comprising  the 
ambient  noise  field,  and  does  not  follow  solely  from  the  inde¬ 
pendence  of  each  of  these  random  sources.9, 14  The  matrix  Q 
is  diagonal  in  the  Kuperman-Ingenito  development  (Ref. 
10,  p.  1991 )  because  of  the  assumption  made  in  their  work 
that  the  spatial  coherence  of  the  noise  depends  only  on  sepa¬ 
ration.  The  matrix  UQU  "  plays  the  same  role  as  K,  in  Bag- 
geroer  et  al.  - 

We  shall  adopt  the  following  model  for  the  CSM  R. 

R=P2Us„s!;U" +  crUU" +p2J  .  (5) 

This  represents  a  single  source  s,„  a  correlated  ( modal )  noise 
component  UQU  "( Q  =  crl),  and  a  white  noise  component 
p2I.  More  general  models,  which  include  a  second  (interfer¬ 
ing)  source,  have  also  been  considered,2  but  Eq.  (5)  will  be 
sufficient  for  our  purposes.  It  is  a  relatively  high  level  of  the 
crUU"  term  that  will  be  of  importance  here.  That  Q  be  diag¬ 
onal  is  not  crucial  to  our  analysis;  only  that  Q  have  full  rank 
M. 

The  first  step  in  any  matched-field  approach  to  estimat¬ 
ing  s()  is  to  construct,  for  each  potential  source  vector 
s  =  s(r,z)  (where  r  and  z  are  range  and  depth),  the  vector 
p(r,z)  =  Us(r,z)  that  the  array  would  have  received.  This 
can  be  done  through  the  use  of  complex  computer  models  of 
the  propagation,  such  as  parabolic  approximation  rou¬ 
tines,11  or  by  employing  the  matrix  U  directly  if  it  can  be 
obtained  from  a  (local)  normal-mode  model. 12  The  conven¬ 
tional  estimator  is  essentially  the  matched  filter: 

CONV(r,z)  =  p(r,z)///fp(r,z)  ,  (6) 

and  the  plot  of  CON  V  ( r,z)  over  ranges  and  depths  of  inter¬ 
est  is  referred  to  as  the  ambiguity  surface.  A  detailed  discus¬ 
sion  of  the  conventional  approach  to  matched  field  process¬ 
ing  is  given  in  Ref.  2.  An  alternative  method  involving,  for 
the  normal-mode  case,  the  passage  to  “mode  space"4  '  is  the 
following  modified  conventional  method: 

MCONV(  r,z) 

=  s (r,z)"(U"U)  'U"RU(U"U)  -  's(r.z)  .  (7) 

It  can  be  shown  that  MCONV  is  the  optimal  linear  processor 
for  minimizing  sidelobe  structure  if  the  distribution  of  po¬ 
tential  sources  resembles  a  noise  component  of  the  form 
crUU".  Equivalently,  MCONV  is  the  optimal  linear  proce¬ 
dure  to  maximize  array  gain  against  noise  having  the  form 
crUU". 

As  with  the  plane  wave  case,  when  the  noise  is  not  white 
the  CON  V  method  is  not  optimal  and  the  optimal  procedure 
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then  involves  a  prewhitening.  This  prewhitening  requires 
knowledge  of  the  actual  noise-only  cross-spectral  matrix. 
Since  such  know  ledge  is  usually  not  available  one  might  pre¬ 
fer  to  use  a  method  that  adapts  to  whatever  the  noise  compo¬ 
nent  happens  to  be.  rather  than  relying  on  an  a  priori  model 
of  the  noise:  the  N1L.  method  is  such  a  data-adaptive  method. 


II.  CAPON’S  MAXIMUM  LIKELIHOOD  METHOD 

Capon's  maximum  likelihood  (ML)  method  is  a  non¬ 
linear  estimator  in  that  the  inverse  of  R  is  required  and  that 
its  resolving  capability  increases  with  SNR.  For  moderate 
levels  of  SNR  it  has  better  resolution  than  CONV.  It  is  con¬ 
venient  to  use  because  it  does  not  require  prior  knowledge  of 
the  ambient  noise  (it  adapts  to  whatever  noise  is  present) 
and  can  be  used  with  any  array  geometry.  When  applied  to 
matched  field  processing  it  often  presents  a  less  ambiguous 
ambiguity  surface  than  CONV.  The  derivation  of  MLcan  be 
extended  w  ithout  difficulty  from  the  plane  wave  case  to  that 
of  more  complicated  propagation  in  a  waveguide  and  it  is 
this  extended  ML  that  we  consider  here. 

For  R  as  in  Eq.  (4)  and  p(r,z)  =  Us(r.z).  Capon’s  ML 
method  leads  to  the  maximum  likelihood  estimator 

ML(r.z)  =  l/[p(r.j)"/?  'p(r,z)\ .  (8) 

A  detailed  discussion  of  this  estimator  and  its  application  to 
acoustic  waveguide  array  processing  is  to  be  found  in  Ref.  2. 
Our  treatment  here  is  different  in  that  we  are  concerned  with 
the  effects  of  a  relatively  high  level  of  modal  noise.  We  find  it 
helpful  to  analyze  the  ML  estimator  in  terms  of  the  eigenvec¬ 
tors  and  eigenvalues  of  R.  The  matrix  R  is  Hermitian  posi¬ 
tive  definite  and  so  has  positive  eigenvalues 
•••  >Ay>  0  and  associated  orthogonal  eigenvec¬ 
tors  z, . zv,  normalized  to  have  norm  equal  to  one.  It  fol¬ 

lows  from  the  form  for  R  given  in  Eq.  (4)  that 

T  vf  .  ,  =  •  •  ■  =  A  N  =  p-  and  that,  for  n  =  M  +  1 . N,  z„ 

is  orthogonal  to  the  columns  of  t/.15'1  If  we  now  rewrite  the 
ML  estimator  as 

ML(«)  =  A,,  '|z"p(r,z)|-)  (9) 

we  see  that  the  terms  from  n  —  \l  +  1  to  /V  are  identically 
equal  to  zero  for  all  p(r,z),  so  they  play  no  role  in  the  ML 
estimator  and  the  denominator  of  ML  reduces  to  a  sum  over 

n  =  1 . V/.  When  phase  errors  are  present  the  eigenvectors 

of  R  are  altered  and  the  terms  corresponding  to 

n  -  M  +  1 . V  do  not  vanish  identically,  causing  the  ML 

estimator  to  degrade.  If  we  eliminate  these  last  .V  —  M  terms 
we  obtain  a  new  estimator,  which  we  call  the  “reduced"  ML 
( RML ). 

In  plane  wave  array  processing  the  steering  vectors 
[  analogous  to  our  p(  r,z)  ]  have  constant  norm,  so  no  norma¬ 
lization  of  ML  is  required.  For  normal  mode  (or  more  gen¬ 
eral)  models  p(r.z),;p(r.z)  is  not  constant  anu  tor  that  rea¬ 
son  it  is  common  to  calculate  a  normalized  ML.  In  our  case 
this  would  involve  multiplying  (8)  by  p  ( r.z )  "p  ( r.z )  or. 
equivalently.  replacing  each  p(r.z)  by  p(/'.r)/ 

|  p(  r.z  '  "pi  r.z  i  | '  .  We  feel  that  this  is  not  a  reasonable  pro¬ 
cedure  and  oiler  analytic  argument  in  favor  of  a  different 
normalization. 


For  the  situation  considered  here  the  noise  is  primarily 
modal  so  that,  in  the  absence  of  a  signal,  the  CSM  R  would 
have  the  form  R  =  VU  "  t-  el,  for  some  small  e  >  0.  In  Ap¬ 
pendix  C  we  show  that  as  e  goes  to  0  the  ML  ambiguity 
surface  for  this  R  goes  to  s( r,z) "s{ r.z).  We  choose  therefore 
to  normalize  our  ML  surfaces  by  dividing  by  s(r,z)"s(r,z), 
instead  of  by  p{r,z)"p(r,z).  This  change  in  normalization  is 
important  when  there  is  a  sizeable  modal  noise  component. 
Figures  1  and  2  illustrate  the  ML  estimator  for  R  as  in  (  5). 
normalized  by  dividing  by  p(r,z)"p(r,z)  and  by 
s(r.z)ns(r.z).  respectively. 


III.  THE  “REDUCED”  MAXIMUM  LIKELIHOOD  METHOD 
(RML) 


The  “reduced"  maximum  likelihood  estimator  (  RML) 
agrees  with  Eq.  (9)  except  that  the  sum  is  over  the  first  M 
terms  only: 


R  M  L  ( r.z ) 


zf/p  (r.z) 


(10) 


We  shall  compare  ML  with  RML  in  our  simulations  and 
observe  the  extent  to  which  ML  degrades  when  phase  errors 
are  present.  To  compute  RML  it  is  not  necessary  to  obtain 
explicitly  the  eigenvalues  and  eigenvectors  of  R:  in  fact, 
RML  can  be  obtained  from  the  matrix  used  in  the  MCONV. 

The  passage  to  "s  space”  that  occurs  in  calculating  the 
MCONV  involves  the  transformation  of  R  to 
T=(U"U)  'U"RU(U"U)  '  and  the  transformation  of 
the  “potential  source”  (or  “steering")  vectors  p  =  p(r.z)  to 
s  =  s(r,z)  =  (U"U)~  'U"p(r,z).  Applying  the  usual 
“maximum  likelihood"  approach  of  Capon,  we  are  led  to 
consider  the  estimator  1/s "T  's,  which  as  we  shall  show,  is 
equivalent  to  the  RML. 

To  show  equivalence  of  the  two  estimators  we  begin 
with  the  eigenvalue-eigenvector  decomposition  of  R 

R=  £  A(1z„z".  (11) 

n  1 


l  K i .  1 .  Ambiguity  surface  produced  using  Ml.  estimator  for  Pckens  wave¬ 
guide  environment  The  level  of  white  noise  is  W  dB  tie.,  a  signal  to 
white  noise  ratio  of  10  dB )  The  level  of  modal  noise  is  20  dB  (  re.,  a  signal  to 
modal  noise  ratio  of  20  dB )  The  steering  vectors  art  normalized  so  that 
p(  r.z) "p(  r.z )  l.forallr.r. 


FIG  2  Ambiguity  surface  produced  using  ML  estimator  for  identical  en¬ 
vironmental  and  noise  conditions  as  in  Fig.  1.  The  steering  vectors  are  nor¬ 
malized  so  that  st  r,z)  ,;s(  r.z)  =  I,  for  all  r.z. 


Since  z„  is  orthogonal  to  the  columns  of  U  for  n>M  +  l.the 

z,,  for  n  =  1 . M  must  be  in  the  linear  span  of  these  same 

columns.1'  '  that  is.  z„  =  Ua„  for  some  a, ,  Hence, 

U"RU  =  y  A„U"z„z!,'U 

II  I 

=  f|,(t/"(/)a„a"((/l,[/).  (12) 

•i  i 

Then  with  fV  =  (  U  "(J) 1 =  W"  we  have 

W  1 V "RUW  1  =  £  A„Wa,X!fV" .  (13) 

«  i 

Since  (  Wa„ )"(  Wa, )  -  a ''W"Wak  =  a"U"t/aA  =  z'X  it 
follows  that  {IFa,,.  r.  =  I  is  an  orthogonal  set,  in 

which  case 

(  W  1  U"RUW  ■')  ■  '  =  W(U"RU)-'W 


=  £  A;  'Wa„^W" .  (14) 

n  I 

By  eliminating  the  W  from  the  left  and  the  right  one  obtains 
(U"RU)  '  =  £  (15) 

ii  i 

Since  s  =  s(r,z)  =  (  U"U)  ~  1  U"p(r,z)  we  have, 
s(r.z)"T  's (r.z)  =  p(r,z)"U(U "RU)  ~ : U "p(r,z)  . 

(16) 


Using  Eq.  (15)  in  Eq.  (16)  we  have 


s (r,z)"T  's (r.z)  =  £  T  „  '  p(r,z)"Uana"U"p(r,z) 

n  I 


=  X  'I*"p(^)I2  •  O7) 

at  I 

Hence. 

RML  =  [s"(r,z)  T  's(z.z)]  ',  (18) 

and  the  RML  estimator  has  precisely  the  same  form  as  the 
ML  estimator  (8).  We  note  that  in  order  to  calculate  RML 
we  need  only  invert  the  single  matrix  U"RU  ( 16),  which  is 
M  by  M 


From  the  above  calculations  we  see  that  RML  can  be 
calculated  as 

RML(r,z)  =  [p(r,z)"U(U"RU)  'U"p(r.z)\  (19) 

The  matrix  (  U"RU)  to  be  inverted  can  be  ill-conditioned 
and  it  is  safer  to  orthogonalize  V  first.  We  can  use  "QR 
factorization"' M>  to  write  U  =  VX,  where  V  is  an  .V  by  M 
matrix  with  M  orthogonal  columns  of  unit  length  and  X  is  an 
upper  triangular  matrix;  this  is  essentially  Gram-Schmidt 
orthogonalization  in  matrix  language.  We  then  have 
V"V  =  /.  and  replacing  U  with  VX  in  (19)  leads  to 

RML  (r.z)  =  [p(r.z)"V(V"RV)  '  V"p(r,z)  ]  "  1  .  (20) 

The  matrix  V"RV  will  be  better  conditioned  than  U"RU 
generally,  and  so  more  accurately  inverted. 

We  see  from  Eqs.  (8)  and  (9)  that  the  ML  estimator 
will  show  a  peak  at  the  true  value  p  =  p(l  =  Us,,  provided  that 
the  denominator  is  close  to  zero  for  p  =  pu  and  larger  for 
neighboring  values.  The  information  about  the  signal  is  then 
carried  by  the  (near)  nulls  of  the  functions  \z'Jp(r.z)  \-.  for 
n  =  2.....JV.  When  the  denominator  is  perturbed  by  the  addi¬ 
tion  of  small  but  nonzero  quantities  these  (near)  nulls  are 
disturbed,  with  the  result  being  the  loss  of  the  peak  at  the 
true  signal  location.  This  is  what  happens  when  the  terms 
corresponding  to  n  =  M  +  1,...,,V  are  no  longer  identically 
zero,  due  to  the  presence  of  phase  errors.  In  our  simulations 
we  compare  the  performance  of  the  ML  and  the  RML  for 
varying  amounts  of  phase  error.  Before  we  pass  to  the  simu¬ 
lations  it  will  be  instructive  to  examine  the  behavior  of  the 
ML  estimator  in  one  special  case,  in  which  U"U  ---  /. 

With  U"U=I  we  have  A,  =  a2  +  cr  +p", 
A2=  -A„=<r+p:,  where  the  signal  power  is 

a 2  =/?:s„"s0  We  can  then  write  the  ML  estimator  as 

ML  =  (  A  ,  '|zf'pj-  +  (cr  +p2)~'  £  |z"p|- 

V  n  =  2 

+  p->  £  |z"p|2)  '  ■  (21) 

As  the  a2  gets  larger  the  value  of  A  ,  1  gets  smaller,  so  the 
first  term  in  the  denominator  gets  smaller,  making  the 
(near)  nulls  of  the  remaining  terms  more  significant  and  the 
peak  of  the  ML  estimator  more  visible.  For  smaller  values  of 
a2,  however,  the  first  term  contributes  a  sizable  nonzero 
quantity  to  the  denominator,  causing  the  ( near )  nulls  of  the 
remaining  terms  to  become  less  significant  and  the  peak  in 
the  ML  estimator  to  drop;  this  is  the  reason  for  the  depend¬ 
ence  of  resolution  on  SNR. 

Holding  a  and  p  fixed,  let  us  consider  the  effects  of  rais¬ 
ing  a.  As  a  is  changed,  the  eigenvectors  z„,  n  —  1 . M 

change,  while  those  for  n  =  M  +  1  do  not.  Ifcr  is  small 
compared  to  a2+p2,  for  n  =  2,...,M,  the  values  of 
U"P(  r,z)  j’  will  be  nearly  zero  for  (r,z)  corresponding  to  the 
true  source  parameters,  but  as  a  increases,  these  (near)  nulls 
will  weaken  and  become  less  distinguished  from  neighboring 
values.  In  addition,  as  a  increases,  A  ,  '  =  (a2  +  cr+p:)  1 
becomes  more  like  A  2  1  =  ,...,  =  A  M '  =  (cr  +  p2)  ~ in¬ 
creasing  the  (unwanted)  contribution  from  the  first  term  in 
the  denominator  of  Eq.  (21 ).  The  overall  effect  is  to  make 
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the  (  near)  null  in  the  signal  direction  more  shallow,  to  raise 
the  background  noise  level  and  to  increase  sensitiv  ity  to  the 
presence  of  the  third  term  in  the  denominator  of  F.q.  (21 ). 

IV.  RESULTS 

In  this  section  we  investigate,  by  means  of  computa¬ 
tional  simulation  within  an  assumed  shallow-water  acoustic 
environment,  the  behavior  of  both  ML  and  RML  estimators 
as  they  respond  to  vary  ing  amounts  of  phase  error  applied 
randomly  to  the  elements  of  the  hydrophone  array.  The 
geoacoustic  environmental  model  used  here  is  the  two¬ 
layered  liquid  half-space  model  of  Pekeris.  Even  though 
the  simplicity  of  the  Pekeris  model  limits  its  general  applica¬ 
bility  .  it  does  possess  features  that  are  v  ery  similar  to  at  least 
one  shallow -water  environment  in  which  matched-field  ex¬ 
periments  have  been  performed. Is  In  this  case,  the  water 
sound-speed  profile  was  almost  isospeed,  and  a  thick  sandy 
sediment  layer  was  present,  which  carried  no  shear  waves 
and  therefore  behaved  like  a  liquid.  The  Pekeris  model  is  also 
a  widely  used  and  well-understood  standard  model.  It 
should  describe  most  of  the  acoustic  characteristics  of  gen¬ 
eral  shallow-water  waveguides.  Since  the  acoustic  wavefune- 
tions  w  ithin  the  waveguide  are  calculated  analytically,  with¬ 
out  need  for  recourse  to  numerical  techniques,  the 
calculation  of  the  pressure  field  vectors  p,  the  transforma¬ 
tion  matrices  L  and  the  modal  amplitude  vectors  s  is  greatly 
facilitated.  This  makes  the  transformation  from  R  to  T.  and 
therefore  the  calculation  of  the  RML  estimator  ( 18).  a  rela- 
tiv ely  simple  matter. 

The  Pekeris  model  consists  of  a  shallow  isospeed  water 
layer  of  uniform  density  overlying  a  faster,  isospeed,  semi- 
infinite  fluid  bottom  of  uniform  (and  usually  higher)  den¬ 
sity  For  a  full  description  of  the  general  Pekeris  model,  the 
reader  is  referred  to  a  standard  text  (c.g..  Ref.  17).  In  this 
work  the  water  depth  is  100  m.  the  sound  speed  in  the  water 
and  sediment  is  1500  and  1621.6  m/s.  respectively,  and  the 
sediment/water  density  ratio  is  1.772.  The  acoustic  source 
has  a  frequency  of  150  Hz.  which  leads  to  the  establishment 
of  eight  propagating  modes  in  this  waveguide.  The  source 
was  placed  at  a  range  of  4000  m  from  the  array  at  a  depth  of 
10  m.  A  vertical  hydrophone  array  consisting  of  21  evenly 
spaced  hydrophones  is  placed  to  sample  the  central  50  m  of 
the  water  column. 

The  difference  in  order  of  the  modal  and  hydrophone 
space  cross-spectral  matrices  here  does  not  lead  to  any  un¬ 
due  computational  problems  in  transforming  from  one  pro¬ 
cessing  space  to  the  other.  Since  the  array  spans  only  50r7  of 
the  water  column  the  matrix  L  "U  ?-/.  and  its  eigenvalues 
will  typically  not  have  the  distribution  described  by  (21). 
However,  many  of  the  broad  features  of  the  eigenvalue  dis¬ 
tribution  will  he  similar  to  t  2 1  ).  which  may  frequently  serve 
as  a  basis  for  interpreting  the  phenomena  observed  w  hen  the 
noise  levels  are  varied. 

White  noise  is  introduced  into  the  system  by  adding  an 
appropriate  constant  |i.e.,/r  in  Eq.  (5)1  to  every  diagonal 
element  of  the  cross-spectral  matrix  in  hydrophone  space  R. 
The  value  of  p  is  determined  by  scaling  against  the  mean 
power  detected  on  the  hydrophones  (equivalent  to  the  mean 
diagonal  element  of  R).  Modal  noise  is  introduced  in  a  simi- 
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lar  way,  by  adding  a  constant  cr  to  the  diagonal  elements  of 
the  modal  space  cross-spectral  matrix  7.  The  value  of  cr  is 
determined  by  scaling  the  diagonal  elements  of  irUU" 
against  the  mean  power  detected  on  the  hydrophones.  In 
these  simulations  the  acoustic  source  level  is  held  steady  so 
that  the  field  detected  at  the  hydrophone  array  due  to  the 
source  remains  constant.  The  values  of p~  and  cr  are  then 
adjusted  to  achieve  the  desired  levels  of  white  noise  and  mod¬ 
al  noise,  respectively.  We  w  ill  quantify  the  values  of/r  and  cr 
in  terms  of  their  dB  level  above  (  +  )  or  below  (  —  )  the 
signal  level.  Note  that  this  method  of  quantifying  the  noise 
levels  ( which  is  convenient  and  helpful  here )  is  the  opposite 
of  that  usually  adopted.  Thus,  a  modal  noise  level  of  —  10 
dB  would  correspond  to  a  signal  to  modal  noise  ratio  of 
-I-  10  dB.  The  same  applies  to  the  relationship  between  the 
level  of  white  noise  and  the  signal  to  white  noise  ratio. 

In  order  to  provide  a  means  of  quantifying  the  rough¬ 
ness  of  the  ambiguity  surfaces,  and  a  measure  of  peak  resolu¬ 
tion.  the  peak  value  ( P ).  mean  background  level  (p).  and 
standard  deviation  of  the  background  ( v)  are  calculated  for 
each  surface.  The  mean  background  level  is  calculated  by 
excluding  a  small  interval  around  the  peak.  Two  measures  of 
peak-to-backgroifnd  resolution  are  used.1''’"  These  are 
PBR1  =  {P  —  p)/p  and  PBR2  =  [P  —  fi)/v.  We  do  not  in¬ 
fer  statistical  behavior  from  these  quantities,  but  merely  use 
them  as  indicators  of  “goodness"  of  the  surfaces. 

Phase  errors  are  introduced  into  the  system  in  the  fol¬ 
lowing  manner.  First  the  maximum  phase  error  is  specified, 
say  0  degrees.  ( We  refer  to  this  number  simply  as  the  “ran¬ 
dom  phase  error"  or  RPE. )  Then  a  diagonal  matrix  D  is 
formed: 

D  =  diag[exp(/<£„ )  ]  .  /t  =  1 . A\  (22) 

where  the  <b„  are  independent  random  variables  uniformly 
distributed  between  —  0  and  +  9.  Then  the  CSM  R  is  re¬ 
placed  by  R  =  DR  D 

Before  we  begin  to  investigate  the  stability  of  the  ML 
and  RML  estimators  to  phase  error,  we  will  reconsider  brief¬ 
ly  the  question  of  estimator  normalization  first  raised  in  Sec. 
II.  The  importance  of  this  question  may  be  illustrated  very 
clearly  if  we  consider  a  case  where  white  noise  is  introduced 
at  a  low  level  compared  with  the  signal  (p~  =  —  30  dB). 
w  hile  the  level  of  modal  noise  is  considerably  higher  than  the 
signal  [cr  =  20  dB).  No  phase  errors  are  yet  introduced.  In 
Fig.  1  we  sec  the  ambiguity  surface  obtained  by  using  the 
summation  expression  (6)  for  the  ML  estimator  with  these 
inputs.  In  this  example  the  steering  vectors  p(r.r)  are  nor¬ 
malized  to  have  unity  norm  [i.e..  p ( r,z )  "p ( r,z )  =  I.  for  all 
r,z\.  Figure  1  indicates  that,  for  the  noise  level  inputs  de¬ 
tailed  above,  this  normalization  procedure  applied  to  the 
ML  estimator  results  in  a  highly  erratic  ambiguity  surface.  A 
ridge  can  be  seen  close  to  the  input  source  position  (10-m 
depth  and  4000-m  range),  but  it  is  masked  by  numerous 
sidelobes.  Several  of  these  have  greater  intensity  than  the 
signal  peak.  Clearly,  reliable  and  unambiguous  estimation  is 
not  possible  on  this  surface. 

Rather  than  normalizing  the  p (r,z)  to  have  unity  norm, 
we  normalize  so  that  their  corresponding  modal  vectors 
s(r,z)  have  unity  norm  [i.e.,  s(r,z)"s(r,z)  =  I,  for  all  r,z ]. 
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With  this  choice  of  normalization  the  ML  estimator  would 
give  a  flat  ambiguity  surface,  in  the  noise-only  case,  if  the 
CSM  is  R  =  LI ' If  we  again  use  the  summation  expression 
(9)  to  recalculate  the  ML  ambiguity  function  for  the 
/r  --  30  dB.  tr  —  20  dB  case  discussed  in  the  previous 

paragraph,  we  obtain  the  surface  in  Fig.  2.  It  must  be  stressed 
here  that  Figs.  1  and  2  have  been  calculated  using  an  identi¬ 
cal  algorithm  [i.e.,  Eq.  (*->)}  with  identical  inputs.  The  only 
difference  between  them  is  the  way  the  steering  vectors 
p(  r.z)  have  been  normalized.  In  fact,  since  the  vector  s (r.z) 
used  in  the  calculation  of  the  RML  estimator  ( 18 )  is  always 
normalized  so  that  s(  r,z)Hs(r.z)  =  1  (for  all  r.z),  and  since 
the  summation  for  the  RML  estimator  ( 10)  is  equivalent  to 
the  summation  for  the  ML  estimator  (9)  (when  the  terms 

for  n  =  3/  -v-  I, . Yare  zero),  we  have  the  surprising  result 

that  the  ML  estimator  (9)  and  the  RML  estimator  ( 10)  are 
equivalent  except  for  the  way  the  steering  vectors  p( r.z)  are 
normalized  in  the  two  cases.  However,  this  is  not  true  in  the 
presence  of  phase  errors,  or  any  other  factors,  which  cause 

the  terms  for  n  —  M  -r  1 . Yin  the  ML  summation  (9)  to 

become  nonzero,  as  we  shall  see  later. 

If  we  inspect  Fig.  2  wc  see  that  the  ambiguity  surface 
contains  an  unambiguous  primary  peak  at  the  input  source 
location  i  10-m  depth  and  4000-m  range):  but  it  is  most  nota¬ 
ble  for  ivvo  other  features.  The  first  is  the  lack  of  contrast 
between  peaks  and  valleys  on  the  surface  itself;  and  the  sec¬ 
ond  is  the  hi«’h  overall  level  of  the  surface,  which  is  com¬ 
pressed  into  the  top  HIT  of  the  available  range  of  estimator 
values.  The  icuson  for  this  may  be  seen  by  inspection  from 
F.qs.  ( 9  >  and  t  21  ).  The  eigenvalues/,  appearing  in  the  sum¬ 
mation  on  the  rhs  of  Eq.  (9)  take  the  values 

/  l=rr  -  tr  -  fr:  A , /  „  ~  ir  +  p:;  /„  .  , /.  N  ~p:\ 

and  since  rr  |>cr '  p:  all  of  these  eigenv  alues  are  dominated 

by  it  and  those  for  n  I . V/  have  similar  magnitude. 

Now.  as  described  earlier,  information  about  the  signal  in 
ML  processing  is  carried  by  the  ( near )  nulls  of  the  functions 

z"p  (r.z)  •’  (for//  --  2 . \f)  in  the  denominator  of  Eq.  (9) 

(the  terms  from //  -  M  -  1  to  A' are  identically  equal  to  zero 
forall  p(  r.z).  so  they  play  no  role  in  the  ML  estimator).  Now 
each  of  the  nonzero  terms  is  multiplied  by  a  factor  ~  1/cr. 
When  tr  is  large  (as  here)  this  factor  tends  to  make  them 
small  over  the  whole  range  of  p(  r.z)  rather  than  at  just  the 
source  position  and.  because  of  the  similarity  of  the  eigenval¬ 
ues  for  n  =  I _ V/  of  comparable  magnitude  to  the  source 

term/.  ,  1  zf'p  [r.z)  This  decreases  the  reciprocal  effect  of 
the  nulls  at  the  source  location  and  leads  ,o  the  shallow  relief 
of  the  details  on  the  surface  and  its  high  overall  background 
level.  Although  the  surface  in  Fig  2  is  very  flat  and  has  a 
high  background  due  to  the  large  amount  of  modal  noise,  it 
does  unambiguously  indicate  the  presence  and  position  of 
the  source  and  is  much  more  useful  than  Fig.  1.  The  reason 
for  this  is  that  it  has  been  produced  using  a  normalization 
procedure  for  the  steering  vectors  which  is  optimized  for  the 
presence  of  modal  noise.  Since  we  w  ill  be  considering  cases, 
in  our  later  simulations,  where  modal  noise  predominates 
over  white  noise,  we  will  use  the  s(  r.z)" s( r.z)  normalization 
throughout. 

We  will  now  proceed  to  consider  the  effects  of  phase 
error.  We  will  be  concerned  here  with  a  case  where  white 


noise  is  introduced  at  a  low  level  compared  to  the  signal 
(p~  =  —  30  dB),  while  modal  noise  is  considerably  higher 
but  still  lower  than  the  signal  (cr  =  —  10  dB).  In  Fig.  3  we 
see  the  ambiguity  surface  obtained  using  the  ML  estimator 
(9)  for  these  input  noise  levels  with  zero  RPE  on  the  hydro¬ 
phones.  This  surface  contains  an  unambiguous  primary  peak 
at  the  true  source  location  ( 10-m  depth  and  4000-m  range). 
The  peak  is  extremely  well  resolved  against  all  the  sidelobes 
present  which  are  minimal  in  intensity.  As  mentioned  ear¬ 
lier.  the  terms  for  n  .V/  +  1  in  thesummation  (9)  should  be 
identically  zero  for  all  p (r.z)  when  no  phase  errors  are  pres¬ 
ent.  so  that  in  this  case  ML  and  RML  are  equivalent  to  each 
other.  Explicit  calculation  of  the  RML  (17)  for  the  same 
input  conditions  confirms  that  this  is  indeed  the  case.  The 
RML  produces  an  ambiguity  surface  identical  to  tuat  of  Fig. 
3. 

The  reason  for  the  high  peak  resolution  on  the  surface  in 
Fig.  3  may  be  seen  by  inspection  of  Eqs.  (9)  and  (21 ).  The 
eigenvalues/,,  appearing  in  the  summation  on  the  rhs  of  ( 9 ) 

take  the  values  A  ,  a«:  4-  tr  +  /r;  A, . /  „  -=■ cr  /r:  and 

/  i;  i . /  \  =  p  .  Since  a:.  <r$>p:  we  may  further  approxi¬ 
mate  these  as/. ,  —  a:  -f  cr;  /.: . A  „  sr  cr; A  ,,  , . /.  x  =  p: 

Passing  to  Eq.  (21 )  we  recall  that  in  ML  information  about 
the  signal  is  carried  by  the  nulls  of  the  functions  z "p(  r.z )  :. 
for  n  ~  2 . V.  In  our  case  here  these  functions  are  identical¬ 

ly  zero  for  n  M  +  1  for  all  p(  r.z)  so  we  are  really  interested 

only  in  the  nulls  of  the  terms  for  n  =  2 . M.  In  (21 )  w  e  see 

that  tr  > cr"  and  therefore  /  ,  1  <A  „  1  (/;  =  2 . V/).  mak¬ 

ing  the  first  term  in  the  denominator  small  and  the  nulls  of 
the  remaining  terms  more  significant.  This  'ead*  to  the  pro¬ 
nounced  signal  peak  observed  on  the  surface. 

Let  us  now  consider  the  effect  of  introducing  phase  er¬ 
rors  onto  the  hydrophones.  Figure  4  shows  the  ambiguity 
surface  produced  using  ML  when  the  RPE  is  30  deg.  We  see 
immediately  that  a  lot  of  degradation  has  taken  place.  Al¬ 
though  the  signal  peak  may  still  be  observed  at  this  point,  it  is 
now  dominated  by  the  sidelobe  structure  and  may  no  longer 
be  unambiguously  resolved.  The  reason  for  this  loss  of  per¬ 
formance  is  primarily  due  to  the  inclusion  of  the  terms  for 


FIG.  3  Ambiguity  surface  produced  using  ML  estimator  for  Pckeris  wave¬ 
guide  environment.  The  level  of  white  noise  is  -  30  dB  (i.e..  a  signal  to 
while  noise  ratio  of  30  dB)  The  level  of  modal  noise  is  -  10  dB  (i.e..  a 
signal  to  modal  noise  ratio  of  10  dB) .  The  RPE  is  0  deg 
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M(j  4  Ambiuuit>  Mirtucc  produced  using  Ml.  eshniaior  for  identical  en¬ 
vironmental  and  noise  conditions  as  in  Fig.  I.  The  RPK  is  30  deg. 

n  Xf  -  1  in  the  summation  for  the  ML  estimator  (9).  With 
no  phase  errors  these  terms  are  identically  zero,  as  we  have 
seen.  However,  the  introduction  of  phase  errors  into  the 
CSM  R.  while  leaving  the  eigenvalues  „  unchanged,  gives 
rise  to  arbitrary  rotation  of  the  eigenvectors  z„  in  the  vector 
space  of  R.  The  projections  of  the  steering  vectors  p ( r.z ) 
along  the  eigenvectors  z„  (//  Xf  +  1 )  then  become  nonzero. 
Although  these  projections  may  be  relatively  small  their  ef¬ 
fect  upon  the  summation  in  (9)  is  amplified  by  the  A  „  1 
factor  by  which  they  are  multiplied.  We  remember  from 
(21 )  that  for  n  Xf  i-  1.  A,,  s  \/p:.  Now  p  is  the  level  of 
white  noise,  which  in  our  case  is  low  (  —  30  dB).  Therefore 
I  //)■’  is  a  large  coefficient  that  can  greatly  exaggerate  the 
contribution  of  the  n  Xf  +  1  terms  to  the  summation,  even 
when  the  RPE  is  comparatively  small.  Using  the  same  envir¬ 
onmental  and  noise  conditions,  we  plot  in  Fig.  5  the  percen¬ 
tage  contribution  to  the  summation  in  the  denominator  of 
the  MI.  expression  (9),  when  rand  z  are  correct,  due  to  the 
ii  Xf  ~  1  terms  We  see  that  the  contribution  rises  rapidly  as 


PHASE  ERROR  (deg) 


FIG  5  Percentage  contribution  of  n  M  +  ]  terms  to  summation  in  the 
denominator  of  expression  f9)  for  the  MI.  estimator.  The  environmental 
and  muse  conditions  are  those  of  the  previous  examples. 


the  phase  error  increases.  It  reaches  50%  with  an  RPE  of  35 
deg  and  95%-  with  an  RPE  of  only  80  deg.  It  asymptotically 
approaches  100%  thereafter.  With  this  in  mind  it  is  easy  to 
understand  how  the  sidelobe  structure  of  the  surface  is  en¬ 
hanced  by  the  presence  of  the  extra,  nonzero  terms  in  the 
ML  summation  (9),  and  the  stability  of  the  estimator  corre¬ 
spondingly  decreased. 

The  RML  estimator  overcomes  this  problem  of  instabil¬ 
ity  due  to  phase  error  by  omitting  the  terms  for  n>A/  -+  1 
from  the  summation.  Some  errors  will  still  be  present  due  to 
the  changes  in  the  projections  of  the  steering  vectors  p(r,z) 
along  the  eigenvectors  z„  for  n^M,  as  these  are  also  rotated 
in  the  vector  space  due  to  the  phase  errors.  However,  this 
effect  is  much  smaller  and  allows  the  RML  estimator  to  re¬ 
main  stable  to  considerably  higher  degrees  of  phase  error. 
Figure  6  shows  the  ambiguity  surface  obtained  using  the 
RML  estimator  with  an  RPE  of  30  deg  (cf.  Fig.  4).  Unlike 
the  ML  case,  the  signal  peak  is  still  unambiguously  localized 
and  sharply  resolved  against  the  sidelobe  background.  For 
this  particular  set  of  environmental  and  noise  conditions,  the 
RML  estimator  remains  stable  up  to  an  RPE  of  85  deg,  when 
the  sidelobes  finally  overcome  the  signal  peak  and  the  esti¬ 
mator  fails.  This  represents  a  remarkable  improvement  over 
ML.  Reference  back  to  Fig.  5  shows  that  for  an  RPE  of  85 
deg  the  n>M  +  1  terms  contribute  more  than  95%  to  the 
ML  summation  (9).  This  is  all  simply  avoided  by  using  the 
RML. 

Figure  7  displays  the  variation  in  PBR1  and  PBR2  as  a 
function  of  RPE  for  the  ML  and  RML  estimators.  The  en¬ 
vironmental  and  noise  conditions  are  unchanged  from  the 
previous  examples.  If  we  individually  compare  the  PBR1 
and  PBR2  curves  for  the  two  cases  we  see  that,  by  both  mea¬ 
sures,  the  RML  estimator  consistently  outperforms  the  ML 
estimator  throughout  the  range  of  RPE  values  investigated, 
except  by  a  small  amount  within  an  interval  between  0  and 
15  deg.  As  the  RPE  is  increased  the  ML  estimator  deterio¬ 
rates  very  rapidly.  The  ML  curves  terminate  at  an  RPE  of  30 
deg,  beyond  which  the  signal  peak  cannot  be  resolved.  The 
RML  curves  also  indicate  a  deterioration  in  performance, 
but  at  a  much  more  moderate  rate.  As  discussed  above,  the 
RML  continues  to  provide  accurate  and  robust  localization 
estimates  up  to  an  RPE  of  85  deg. 


FIG.  6.  Ambiguity  surface  produced  using  RML  estimator  for  identical 
environmental  and  noise  conditions  as  in  Figs.  I  and  2,  The  RPE  is  JO  deg. 


2499 


J  Acoust  Soc  Am..  Vol  87,  No  6,  June  1990 


Byrne  et  al. :  Data-adaptive  matched-field  processing 


2499 


PHASE  ERROR  (deg) 

F  IG.  "V  Variation  ot'PBRl  and  PBR2  as  functions  of  R PE  tor  both  ML  and 
'  Ml.  estimators.  The  emironmentul  and  noise  conditions  are  those  of  the 
previous  examples.  The  four  curves  are:  (a)  RML(PBR2) — shaded  dia¬ 
monds:  l  b)  RML(  PBR1) —  unshaded  diamonds;  te)  MLtPBR2)  — 
shaded  squares;  Id)  ML  I  PBR1 ) — unshaded  squares. 


V.  CONCLUSIONS 

A  modification  of  the  data  adaptive  maximum  likeli¬ 
hood  (ML)  method  is  presented  that  shows  substantially 
improved  stability  in  matched  field  processing  for  source 
localization.  Instabilities  in  the  ML  estimator  arise  from 
phase  errors  in  the  signals  received  at  the  hydrophone  array. 
These  may  be  caused  by  a  number  of  physical  and  environ¬ 
mental  factors.  We  simulate  them  here  by  introducing  ran¬ 
dom  phase  errors  at  the  hydrophones.  Applying  the  eigen¬ 
vector/eigenvalue  decomposition  of  the  cross-spectral 
matrix  R.  the  ML  estimator  is  expanded  as  the  reciprocal  of 
a  sum  of  terms,  one  corresponding  to  each  eigenvector.  In 
theory,  most  of  these  terms  should  be  zero  when  the  source 
has  been  correctly  located,  leading  to  a  peak  in  the  estimator. 
Due  to  perturbations  in  R  caused  by  the  phase  errors,  how¬ 
ever,  terms  that  should  be  zero  across  the  entire  ambiguity- 
surface  are  not.  leading  to  a  drop  in  peak  level.  By  removing 
these  erroneous  terms  the  estimator  may  be  stabilized.  We 
term  the  new  estimator  that  results  from  this  procedure  the 
“reduced  maximum  likelihood"  or  RML  estimator.  Dimen¬ 
sionality  reduction  procedures  are  used  to  calculate  the 
RML.  avoiding  inversion  of  the  full-sized  matrix  or  the  cal¬ 
culation  of  eigenvalues. 

Using  a  Pekeris  waveguide  model,  simulations  of  acous¬ 
tic  propagation  from  a  single  source  to  a  21 -hydrophone  re¬ 
ceiving  array  were  used  to  examine  the  RML  and  ML  esti¬ 
mators.  In  our  particular  model  the  waveguide  supports 
eight  modes  and  application  of  the  RML  estimator  reduces 
the  full  problem  of  dimension  21  to  one  of  dimension  8.  Con¬ 
sidering  an  environment  with  a  large  amount  of  modal  noise 
(the  signal  to  modal  noise  ratio  is  —  20  dB )  and  little  white 
noise  (the  signal  to  white  noise  ratio  is  30  dB),  it  is  found 
that  using  an  ML  estimator  which  is  optimized  for  the  pres¬ 
ence  of  modal  noise  ( by  performing  an  appropriate  normali- 
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zation  of  the  steering  vectors)  produces  a  far  superior  and 
more  useful  ambiguity  surface  than  using  an  ML  estimator 
which  is  optimized  for  the  presence  of  white  noise.  In  the 
absence  of  phase  errors  the  ML  and  RML  estimators  are 
equivalent  apart  from  the  steering  vector  normalization 
schemes  they  employ. 

In  the  case  of  an  environment  with  a  significant  amount 
of  modal  noise  ( the  signal  to  modal  noise  ratio  is  1 0  dB )  and 
little  white  noise  (the  signal  to  whitenoise  ratio  is  30dB)  the 
two  methods  are  compared  as  they  respond  to  increasing 
degrees  of  phase  error  introduced  on  the  hydrophones.  The 
results  obtained  clearly  demonstrate  the  superiority  of  the 
RML  estimator  over  the  ML  estimator.  The  ML  estimator 
deteriorates  rapidly  with  phase  error,  due  to  the  contribu¬ 
tions  of  large  (previously  zero)  terms  in  the  summation 
expression  for  ML.  The  ML  estimator  fails  completely  when 
the  random  phase  error  is  greater  than  30  deg.  The  RML 
estimator  avoids  this  problem  by  omitting  the  destabilizing 
terms  from  the  summation.  The  terms  remaining  in  the 
RML  summation  are  also  affected  by  the  introduction  of 
phase  errors,  but  to  a  much  lesser  extent.  The  RML  estima¬ 
tor  continues  to  perform  accurately  and  reliably  up  to  a  ran¬ 
dom  phase  error  of  85  deg. 
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APPENDIX  A:  FORMULATION  OF  p  =  t/s:  LOCAL 
NORMAL-MODE  MODEL 

In  Sec.  I  we  discussed  the  normal-mode  model  for  prop¬ 
agation  and  showed  that  the  single  snapshot  ( narrow-band ) 
data  vector  corresponding  to  a  single  source  had  the  form 
p  =  Us.  where  the  N  XM  matrix  U  is  independent  of  the 
source  and  depends  only  on  the  array  and  the  acoustic  envi¬ 
ronment.  All  source  dependencies  are  in  the  vector  s,  whose 
entries  are  given  in  Eq.  (3).  However,  the  typical  normal¬ 
mode  approach  has  the  deficiency  of  being  valid  only  for 
depth-dependent  media,  and  does  not  permit  inclusion  of 
range  dependencies.  In  this  Appendix  we  shall  show  that  the 
p  —  formulation  and  preceding  results  are  valid  for  com¬ 
plicated  environments  where  the  ocean  channel  has  sound 
speed  and  density  exhibiting  both  range  and  depth  depend¬ 
encies,  as  well  as  a  rough  surface.  Such  assumptions  in  mod¬ 
eling  an  ocean  channel  can  be  included  in  an  extension  of  the 
normal-mode  method  to  obtain  solutions  in  terms  of  “local" 
normal  modes.  We  shall  use  well-known  methods"  :i  in  de¬ 
riving  solutions.  First  we  incorporate  density  changes  into 
effective  sound  speed  profiles,  and  then  we  obtain  solutions 
in  terms  of  orthogonal  eigenfunctions  and  coupled  modal 
amplitudes. 

Consider  a  sound  source  positioned  below  the  origin  of  a 
cylindrical  coordinate  system  on  the  z  axis  at  z  =  z, .  Assu  - 
ing  a  cylindrically  symmetric  ocean  channel,  the  surface  is  at 
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Z  =  d(r).  and  the  bottom  of  the  ocean  sediments  is  at  2  =  D. 
The  ocean  is  assumed  to  have  range-  and  depth-dependent 
sound  speed  <.'(/•, z)  and  density  p„(r.z).  Let  the  acoustic 
pressure  be  p(r.z.i).  Assuming  harmonic  time  dependence 
with  frequency  to  radians  per  second,  we  have 

p{r.z.n  =  P(r.z)c  .  ( A 1  > 

We  suppress  any  time  dependence  in  the  surface  [choosing 
to  view  d(r)  as  a  member  of  an  ensemble  of  random  func¬ 
tions  remaining  fixed  in  source-receiver  travel  times]  and 
hence  also  in  P.  which  will  satisfy  the  differential  equation 

'V/3)  +  A  :P  =  [  —  S(r)’S(z  —  2, )  ]/2rrr .  (A  2) 

along  with  the  horizontal  boundary  conditions, 

P{r.d(r))  —  0  (A3) 


and  where  s  is  a  source  dependent  vector  with  the  M compo¬ 
nents 

.v,„  =  <!»,„ (r)  ,  /m  =  1.2 . M.  (A14) 

Note  that  the  ,v„,  depend  on  source  range  and  depth.  In  the 
case  when  there  is  modal  and  white  noise,  in  addition  to  a 
particular  signal  of  interest,  one  writes 

x  =  f/,  +  t/y  +  T]  ,  ( A 1 5 ) 

as  in  Eq.  (1).  Even  though  it  appears  that  Eq.  (A13)  de¬ 
pends  on  source  range  ( since  that  is  where  the  origin  is  posi¬ 
tioned  )  it  is  actually  not  dependent  on  the  distance  between 
the  source  and  receiver  but  on  the  actual  location  of  the 
receiver;  i.e.,  for  different  source  locations  the  U,„„  are  still 
the  same. 


dP 

—  (/\£>)  =0.  ( A4) 

dz 

Equation  (A3)  corresponds  to  a  pressure  release  surface, 
while  Eq.  ( A4)  is  descriptive  of  a  perfectly  hard  subbottom. 
In  Eq.  ( A2),  A  =  to/c  is  the  local  wavenumber.  By  defin- 


P(r.z)  =  \pt)(r,z)  P'(rsz)  ,  (A5) 

and  using  Eq.  (A2),  P‘  satisfies 

V~P'  +  k2.P'  =  [  —  8(r)-S(z  —  2, )  ]/2 nrp,, .  ( A6) 

where 

ki  =k2  +  lp!rV-(Pl)  ,  :Vptl)  (A7) 

is  an  effective  local  wavenumber  Hence,  by  defining  an  “ef¬ 
fective”  sound  speed,  c,  =  <o/k,.,  one  may  implement  den¬ 
sity  variations  in  numerical  routines  via  effective  sound 
speed  profiles  and  modified  source  strength. 

To  derive  coupled  modes  we  let"  :i 

P'(r.z)  =  <t>„(r)^„(z;r)  ,  (A8) 

n  I 

where  (z\r)  are  solutions  to  the  local  normal-mode  prob¬ 
lems 


-r-  +  [A;’(r,z)  -*;(/•)]  ¥„  =  0 , 


>1 >„(d(r)\r)  =  0 


-(p,  ’*„) .-  „  =  0. 


(A10) 


(All) 


See  Ref.  11  for  the  equations  that  uniquely  determine  the  <t>„. 
Now,  if  we  use  the  typical  definition  for  the  vector 

p=  [P(r.z,  ),P(r,z2),P(r,zt ),...,  P(r,zs  )]'  ,  (A12) 

and  assume  only  a  finite  number  (M)  of  modes  propagate, 
then  p  =  (/, ,  where  f/ is  an  N  XM  matrix  with  entries 

n  =  1,2 . N, 

U =  xp„( r.z„  )  T*,,,  (z„:r)  ,  ( A 1 3 ) 

m  =  1,2  ,...,M, 


APPENDIX  B:  FORMULATION  OF  p  =  US:  PARABOLIC 
APPROXIMATION  MODEL 

For  actual  applications,  the  system  of  equations  in  the 
preceding  appendix  may  be  quite  difficult  to  solve.  The  Para¬ 
bolic  Approximation  method  simplifies  the  modal  ampli¬ 
tude  equations.  If,  instead  of  using  Eq.  (A5),  one  uses  the 
assumption"  that 

P(r,z)  =  ^p„(r,z)  u(rj)H'n"{knr)  ,  (Bl) 

where  in  the  farfield 

k(,r)  ~  Jl/iirkjr  e'Kr ,  A,/>  1  ,  ( B2 ) 

then,  away  from  the  source,  u{r,z )  satisfies  the  parabolic 
approximation 

-^-7-  +  2 ik0^-+  [A  J(r,z)  —  A p,  ]  u  =  0  ,  (B3) 

dz-  dr 

where  A„  =  (o/c,}  is  a  reference  wavenumber  and  c()  is  a  refer¬ 
ence  sound  speed.  Now  if  one  writes 

u(r,z)  =  £  r„(r)'P„(z;r)  ,  (B4) 

n  ----  I 

with  'P,,  determined  as  in  Eqs.  ( A9)-(  A1 1 ),  then  T„  will  be 
solutions  of 

dr 

— - rr- [«'»('■)  +  Ao]E„=-  £.4 . (OT,,  , 

dr  2A0  1 

(B5) 

where  A„,„  are  integral  functions  of'P,,  and  'P,„,  given  in  Ref. 
11. 

If  a  is  a  vector  of  length  M  whose  components  are  modal 
amplitudes  T„,  and  if  we  approximate  the  solution  by  as¬ 
suming  M  modes  propagate,  then  a  satisfies  the  first-order 
problem 

a'  =  ll(r)a,  (B6) 

where  11  is  an  M  X  M  matrix  with  entries 

rf-  (*7.  +  A  f, )  -A . (r) ,  m  =  n, 

11 =  I  2A„  (B7) 


-A . (r) 


m^n. 


Hence,  if  exp[  J,r,ll(r)  dr]  is  the  matrix  exponential  corre¬ 
sponding  to  the  integral  of  11  then  one  can  write 
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a  =  exp^  I  II(r)drjc.  (B8) 

where  c  is  a  constant  vector  determined  by  the  initial  condi¬ 
tion. 

The  advantage  to  this  method  is  that  determination  of 
the  matrix  exponential  in  Eq.  ( BS )  is  perhaps  easier  than 
solving  the  system  of  coupled  equations  in  the  previous  for¬ 
mulation.  In  this  ease,  the  p  —  L\  is  still  valid,  with 

s  =  a// (k[tr)  .  (B9) 

It  is  important  to  again  note  that  (’  is  only  dependent  upon 
receiver  locations  and  the  particular  ocean  characteristics  at 
the  receiving  array,  while  s  contains  all  source  information. 

APPENDIX  C:  ESTABLISHING  THE  ML  ESTIMATE 
WHEN  R—UUH 

To  establish  the  existence  of  the  limit  of  the  denomina¬ 
tor  of  the  maximum  likelihood  estimator  [i.e.. 
p(r.z)"(  UU  "  4-  el)  '  p  ( r.z )  ]  as  e  —  0  we  first  write 

(UU"  +  eI)U=  U(U"U+eI)  .  (Cl) 

from  which  it  follows  that 

(  UU  "  -  el)  '  U  =  U(  U  "U  +(!)'.  (C2) 

Now  we  write  the  denominator  of  the  ML  estimator  as 
p (r.z)"(  UU  "  +  el)  1  p (  r.z)  and  using  the  relation 
p  (r.z)  -—  Us(r.z)  we  have 

f>(r.z)"(UU"  +  el)  1  p  ( r,z ) 

=  s (r.z)"U"(UU"  +  el)  '  'Us(r.z)  ,  (C3) 

and  using  (Cl )  we  find 

p(r,z)"(  UU"  +  el)  'p(r.z) 

=  s(r.z)"U"U(U"U  +  el)  's  (r.z) 

=  s(r.z)"(U"U  +  el)(U"U  +  el)  's  (r.z) 

-  es(r.z)"(  U"U  +  el)  's (r.z) 

=  s(r,z)"s(r.z)  -  fs(r,z)"(  U"U  +  el)  's (r.z)  . 

(C4) 

Since  s(r.z)"(  U"U  +  ef)  's (r.z)  goes  to 


s  (r,z)"(U"U)  's(r.z)  as  e  goes  to  zero,  it  follows  that  the 
entire  expression  goes  to  s(r,z)"s(r,z) . 
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